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AW  INCREMENTAL  CHARGE  METHOD  FOR  THE  ANALYSIS  OF  NONLINEAR  CIRCUITS 

Stephen  John  Nuspl^  M.S. 

Department  of  Electrical  Engineering 

University  of  Illinois,  196^^ 

Conventional  numerical  methods  for  nonlinear  transient  circuit  analysis 
apply  standard  numerical  procedures  to  the  equations  obtained  from  Kirchhoff's 
lawSc   In  this  thesis  is  developed  a  method  of  analysis  in  which  only  one  of 
Kirchhoff's  laws  is  forced  to  hold  and  the  other  is  used  to  make  num.erlcal  cor- 
rections»   The  method  developed  here  assumes  that  the  current  law  holds  and  the 
voltage  law  may  he  violated „ 

To  obtain  sufficiently  simple  numerical  formulas;,  a  circuit  with  a 
tree-type  topology  is  assumed.   Relying  on  this  the  basic  charge  distribution 
and  node  voltage  calculation  formulas  are  derived  for  the  general  branch  which 
may  consist  of  a  resistance,  a  capacitance,  an  inductance,  a  voltage  source  or 
any  series  combination  of  these „   Since  this  method  favors  current  controlled 
nonlinear  elements,  an  analogous  but  more  restricted  method  for  conveniently 
including  voltage -controlled  devices  is  developed,   A  brief  discussion  of  errors 
and  their  compensation  and  of  stability  is  included. 

Because  of  the  importance  of  diodes  and  transistors,  the  technique  by 
which  they  can  be  efficiently  analyzed  by  the  method  is  developed.   The  advan- 
tages are  that  the  method  provides  for  the  models  the  means  for  making  numerical 
corrections e   A  technique  by  which  transmission  lines  can  be  handled  is  also 
described, 

A  program  based  on  this  method  is  presented  through  a  brief  description. 
To  illustrate  the  use  of  the  method  four  sample  circuits  and  the  waveforms  they 
produced  are  given.   It  is  concluded  that  the  method  has  some  decided  advantages 
which  should  prove  to  make  it  useful,  but  not  enough  experience  has  been  had 
[With  it  to  make  any  meaningful  comparison  with  other  efficient  methods  for  non- 
linear  circuit  analysis. 


1,   INTRODUCTION 

In  the  analysis  of  transients  in  nonlinear  circuits  one  usually 
attempts  to  obtain  a  first-order  solution  by  the  use  of  a  linear  transform 
method  which  relies  on  a  piecewise  linear  approximation  of  the  nonlinear 
elements.   When  this  is  no  longer  accurate  enough  or  leads  to  cumbersome 
algebraic  manipulations,  a  more  basic  approach  must  be  used.   This  normally 
consists  in  the  writing  of  the  differential  equations  of  the  circuit  obtained 
from  Kirchhoff 's  laws  and  manipulating  these  until  a  set  of  first-order  linear 
equations  in  the  form 


X, 


kjm,  1   L  ^-'- 


=   a. 


m,  n 


X. 

1 


n,l 


is  obtained.   a   is  either  a  constant  or  a  variable  depending  on  any  x. „   For 

any  but  the  simplest  of  circuits  a  numerical  integration  must  be  used.   This 

method  has  a  number  of  disadvantages  which  become  more  and  more  severe  as  the 

size  of  the  circuit  increases:   finding  the  coefficients  a,  .  is  usually  a  time- 

ki 

consuming  job;  the  evaluation  of  the  coefficients  at  each  step  requires  con- 
siderable computer  time  since  the  coefficients  in  general  contain  some  terms 
related  to  the  nonlinear  elements;  the  addition  of  only  one  more  element  in  the 
circuit  might  require  that  new  expressions  for  all  coefficients  be  found  and 
programmed.   These  inconveniences  severely  restrict  the  practical  use  of  this 
technique  for  the  solution  of  transients  in  larger  nonlinear  networks „ 

Mechanizations  of  the  above  method  which  essentially  separate  the 
■  inear  elements  from  the  nonlinear  are  found  in  the  literature.  ^' ^^   The  linear 
'art  of  the  circuit  is  usually  described  by  a  coefficient  matrix  which  is  auto- 
latically  assembled  and  inverted  by  program.   The  linear  and  nonlinear  portions 
•re  then  integrated  separately  and  matched  at  convenient  time  intervals.   This 
lass  of  methods  is  excellent  in  that  only  a  description  of  the  circuit  topology 
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and  elements  must  be  supplied.   However,  considerable  programming  effort  is 
required  to  make  the  methods  functional  and  sometimes  there  are  stability 

problems . 

The  Incremental  Charge  (IQ)  method  described  in  this  thesis  in  essence 
goes  back  one  step  to  more  basic  principles.   In  all  of  the  above  methods  it 
was  assumed  that  Kirchhoff 's  current  and  voltage  laws  hold  and  from  these  the 
differential  equations  were  obtained.   In  the  Incremental  Charge  method  we 
impose  one  of  Kirchhoff 's  laws,  allow  the  other  to  be  temporarily  relaxed  and 
use  the  deviations  from  this  law  as  a  means  of  making  numerical  corrections. 
The  process  must  always  tend  towards  the  reduction  of  the  deviations.   Programs 
for  the  dc  analysis  of  circuits  ty  the  use  of  these  principles  have  been  success- 
fully used'-^'^-'  but  the  extension  of  these  principles  into  the  time  domain  does 
not  yet  seem  to  have  been  attempted. 


2.   DEVELOPMENT  OF  THE  INCREMENTAL  CHARGE  METHOD 

The  basic  philosophy  behind  the  IQ,  method  is  that  one  of  Kirchhoff 's 
laws  will  be  allowed  to  be  violated  in  order  to  permit  numerical  incremental 
corrections o   We  can  allow  either  the  voltage  around  a  circuit  loop  to  be  in 
error  by  a  small  amount  or  we  can  relax  the  requirement  that  the  algebraic  sum 
of  the  currents  flowing  into  a  node  be  zero„   It  was  found  that  allowing 
Kirchhoff 's  voltage  law  (KVL)  to  be  violated  and  strictly  enforcing  Kirchhoff 's 
current  law  (KCL)  seemed  to  result  in  a  more  convenient  formalization  of  a 
numerical  procedure.   It  is  therefore  this  principle  on  which  the  IQ  method  is 
based, 

2=1  Circuit  Ordering  and  Topological  Restrictions 

The  first  task  is  to  devise  a  convenient  scheme  for  presenting  the 
circuit  topology  to  the  computer  and  an  order  among  the  circuit  elements  which 
the  program  follows  during  the  execution  of  the  method.   To  accomplish  this  each 
aode  of  the  circuit  will  be  assigned  a  number  which  is  also  defined  as  the 
"level"  of  that  node.   This  notion  of  levels  of  nodes  is  introduced  here  only 
for  the  reason  that  it  conveys  better  the  concept" of  ordering  among  the  nodes 
than  the  word  "number,"  The  reference  node  is  always  assigned  the  number  1„ 
Phese  definitions  are  illustrated  in  Fig,  2,1,   Between  any  two  nodes  there  can 
38  branches  which  may  consist  of  a  capacitance,  a  resistance,  an  inductance,  a 
Aoltage  source  or  any  combination  of  these.   These  branches  also  have  numbers 
issigned  to  them  but  since  only  the  knowledge  of  which  ones  are  under  a  partic- 
|ilar  node  is  required,  the  ordering  can  be  arbitrary. 

To  define  a  sufficiently  simple  numerical  procedure  a  number  of 
j-estrictions  will  be  placed  on  the  circuit  topology.   First  it  will  be  assumed 
■hat  the  topology  is  in  the  form  of  a  tree  as  shown  in  Fig,  2,1,   The  forcing 
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Figure  2a,      Assumed  Form  of  Topology  for  Development  of  Method. 


-5- 

function  to  which  the  circuit  must  respond  is  in  the  branch  at  the  apex  of  the 
diagram  and  in  the  form  of  a  voltage  source.   The  second  restriction  is  that 
there  be  at  most  one  branch  between  two  nodes  of  which  the  reference  node  is 
not  one»   These  restrictions  may  seem  drastic  at  present,  but  many  nonlinear 
circuit  problems  fall  directly  into  this  type  of  topology  and  many  others  which 
do  not  conform  to  these  assumptions  can  be  handled  but  the  error  will  be  slightly 
larger  and  other  conveniences  of  the  method  are  forfeited „   The  second  restric- 
tion will  eventually  be  eliminated  entirely  and  the  relaxation  of  the  first  will 
be  discussed  later, 

2.2  The  Charge  Distribution  Criterion 

There  are  essentially  two  steps  in  the  numerical  procedure.   In  the 
first  a  current  is  injected  into  the  highest  node  and  passed  to  successively 
lower  level  nodes  via  the  branches  under  each„   When  the  current  distribution 
is  completed,  the  currents  in  all  branches  are  known,  enabling  us  to  calculate 
the  voltage  across  each  element.   The  second  step  then  consists  in  starting 
from  the  reference  node  and  determining  what  each  successively  higher  level  node 
voltage  will  be. 

To  develop  this  method,  let  us  consider  the  general  node  of  Fig„  2.2 
which  is  assumed  to  have  n  branches  immediately  below  it.   In  a  time  interval 
At^  a  current  I   or,  equivalently,  a  charge 

is  injected  into  node  k  from  higher  level  nodes.   By  a  distribution  process 
which  is  to  be  developed  branch  j  receives  an  amount  of  charge  AQ..,   The 
assumption  that  KCL  holds  implies  that 
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Figure  2„2o      The   General  Node „ 


■7- 


AQ^     =     ZaQ  (2.1) 

The  branches   in  Fig.    2.2  are  not   shown  connected  to  node  k  to  emphasize   that 
voltage  V     and  all  v..  's  may  have   different  values   if  ICVL  does  not   hold        If  we 

K.  J 1 

required  that  it  does  hold,  we  have 

+     +  +  ,    , 

V,  .  =  V-,  .  =  V^  .  =  „  „  „  =  V  (?    p) 

ki    li    21         ni  \^'^J 

since  the  points  (Y)  to  (n^   at  the  branch  tops  are  connected  to  node  k  in  the 
actual  circuits,   If  we  perm.it  this  last  set  of  equations  to  be  violated,  we  are 
allowing  voltage  differences 

Av  . .  =  V,  .  -  V .  .  ( 2  o  ^ ) 

Ji    ki    ji  v-o; 

If  this  Av   approaches  zero,  KVL  will  again  hold.   Hence  an  exact  solution  of 
the  circuit  will  be  had  if  this  condition  holds  over  the  entire  network.   There- 
fore, the  object  of  the  numerical  method  will  be  to  keep  this  voltage  difference 
as  small  as  possible,  or  at  least  below  a  tolerable  level. 

The  problem  now  is  to  distribute  the  charge  AQl.    among  the  branches  in 
such  a  manner  that  all  branch  voltages  v.  will  experience  an  approximately  equal 

J 

increase.   To  accomplish  this  we  define  a  set  of  differences  by 

AQ..  =  Q..  -  Q.,   ^ . 

A^Q  .  =  AQ..  -  AQ.,.   .  (2,^) 

A^Q   =  A^Q..  -  A^Q...  ^. 
Ji      Ji      j(i-l) 

The  voltages  developed  across  the  three  types  of  elements  are  then  given  by 
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Q.. 


cji 


V      =  li 


2 


('A+  '> 


>.2 


(2.5) 


Since  we  are  allowing  each  branch  to  contain  all  three  types  of  elements^  we 
see  ^hat  only  A^Q. .  can  he  changed  arbitrarily;  the  other  charge  differences 

J 

must  be  determined  by  equations  2.k.      Also  define 


A  +     +     + 

Ji    Ji    j(i-l) 


(2.6) 


We  are  now  in  a  position  to  define  three  series  impedance  factors  between  point 
(T)  of  Figo  2,2  and  the  reference  node: 


Av  ^..At. 
cEjx  1 


cj 


AQ. 


(2„7) 


Ji 


Av^„..At. 


Rj     A^Q. 


(2.8) 


Ji 


Av^^.  At. 
LEji   1 


^'j    a\. 


(2.9) 


Ji 


where  v  ,  v   and  v^.^  are  the  voltages  developed  across  these  series  impedance 
cE   RE      -LiE 

factors  and  which  satisfy  the  relations 


^ji  ^  ^cEji  ^  ^REji  +  ^LEji 


(2.10) 


and 


Av   =  Av    .  i-  Av„^  .  .  +  A\'-T  ^, . . 
"ji     cEji     RF.J1     LEji 


(2.11) 


We  can  also  define  a  total  branch  impedance  factor 
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z .  =  z  .  +  z^ .  +  z^ . 

J     cj    Rj     Lj 


(2.12) 


and  two  admittance  factors 


Y  =-^ 


(2a3) 


n 


Y,  =  Z  Y. 


(2aif) 


Phese  factors  serve  as  a  means  of  predicting  branch  voltage  increases  for  a 
5iven  change  in  AQ.  or  in  one  of  the  higher  differences » 

Assume  that  the  numerical  procedure  has  been  under  way  for  a  number  of 
:ycles  and  that  after  cycle  (i-l)  is  completed  there  is  a  residual  voltage  error 
^^j(i-l)  ^^"^^^^^  "°ds  ^  ^^d  branch  j„   The  first  task  is  to  raise  the  branch 
roltage  v   to  the  node  voltage  v^^^,^)  ^y  injecting  an  appropriate  amount  of 
charge  which  we  shall  define  as  A^Q..  into  the  brancho   Assuming  that  the  three 
-mpedance  factors  are  known  for  the  branch,  this  condition  requires  A^Q..  be 
such  that 


AQ 


^(i-D^^^Ki-D^^Si),   ,  ^^'S(i-i 

At.  ^cj  ^         AT 


^Rj    At,  \.j 


Av  . /  .  .  ^ 

j(l-l) 


his  yields 


(2.15) 


^      At. 

J 


Av  -  z      i^^Jiiiill     rz     +  z    )  "^hilill 

L   j(i-l)     cj  \        At.   )-  ^^cj  ^  ^Rj^  ^aE~~. 


(2=16) 


his  is  the  basic  "correction"  formula  in  the  method  and  will  be  discussed  in 
ore  detail  later.   After  determining 
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for  every  branch  and  summing,  we  have  distrihuted  an  amount  of  cb_arge 


Wi  -  .|  Vji 


(2.18) 


This  leaves  an  additional  charge  to  be  distributed: 

^\i  =  ^\i  -  w^  ^'-^^ 

Since  it  is  assumed  that  at  this  point  all  branch  voltages  are  at  the  level 

V  /    X,  this  additional  charge  must  be  distributed  so  that  all  branch  voltages 
k(i-l)' 

will  increase  the  same  amount.   If  branch  j  receives  an  additional  A^Q^^  then 
Since  by  definition 


4Si  =  ''^(i-i)  *  4si 


^«ji  =  ^'^Jd-D  "  4Si 


(2.21) 


we  get  the  result 


4Si  =  4Si  =  Vji  '^■^^' 
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The  usefulness  of  this  lies  in  the  fact  that  the  additional  charge  A  Q   will 

A  ji 

cause  an  increase   in  branch  voltage   of 


Vji        ^cj      At.      ^  ^Rj      At.      ^   hj      At.      -     At.       ^^cj    ^  ^Rj    ^  ^Lj  ^    =  Y;^ 


J      1 
(2.23) 


3r 


A  Q..    =  Y.At.A.v^.  (2.2^) 

/hich  gives  the  amount  of  charge  which  must  be  injected  into  branch  j  to  raise 

Lts  voltage  ^^v   .   We  want  the  voltage  of  all  n  branches  and  that  of  the  node 
:o  increase  the  same  amount : 

"^A^ji  ^  ^A\i  ^  constant  for  j  =  1  to  n  (2„25) 


fence 


herefore 


AQ..    Y. 
A\i    ""k 


^Si=^A\iXY:^  (2.28) 
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This  is  the  desired  forimla  for  distributing  the  charge  ^p9^^   among  the  branches, 
The  term 


I 


\   At.Av 

1  A  Kx 


is  in  the  form  of  an  admittance  and  since  the  expression  contains  only  terms 
related  to  node  k^  we  could  define  it  as  the  "admittance  factor"  for  node  k„   A 
node  impedance  factor  can  therefore  be  defined  by 

2o3  Node  Voltage  Calculation 

The  above  criterion  is  used  to  channel  the  charge  injected  into  the 
top  node  through  all  the  branches  and  other  nodes  and  into  the  reference  node o 
When  this  is  completed,  all  the  branch  currents  are  known  and  therefore  the 
voltages  across  them  can  be  calculated.   The  remaining  step  is  to  calculate  the 

node  voltages  v,  . , 
^    ki 

Through  possible  errors  in  the  impedance  factors  all  the  branch 

voltae-es  v"^  will  not  have  the  same  value „   Therefore  some  sort  of  mean  value 
ji 

for  V   must  be  defined.   There  are  a  number  of  possible  ways  of  determining 
ki 

this  mean,  the  most  primitive  of  which  is  the  average  value : 

V  .  =  i  Z  vt.  (2o30) 

ki   n      ji 

A  second  formula  may  be  obtained  by  imposing  the  conservation  of  energy  conditioi 
at  node  k,  which  implies  that  no  energy  is  lost  through  the  niMierical  process o 
This  condition  states 


-13- 


^i'^S.l  =    .^,   ^Il-^Sl  (2.31) 

J=l 


and  hence  we  have 


n      AQ.  . 

^-^5rji^  (-32) 


This  is  effective  when  AQ   and  all  AQ..'s  have  the  same  sign„   But,  if 

n 

|ao  .  I  «  E  Iaq..  I 

which  occurs  if  one  of  the  branches  under  node  k  contains  a  dc  source,  numerical 
instabilities  may  result.   A  modified  version  of  the  formula  is 

n 

Z  vT.  Iaq.  .  I 

j=i  J"   J"  ,   , 

\^-~ (2.33) 

Z  |aq..  I 

When  all  ^Q . • ' s  have  the  same  sign  energy  is  again  conserved;  otherwise  the 
node  voltage  will  be  the  average  between  the  branches  delivering  power  and  those 
receiving  it. 

In  a  third  method  the  impedance  factors  defined  for  the  current 
distribution  are  used: 

n 

Z  Y.vT. 
1=1  J  J^ 

k 

["his  has  the  effect  of  latching  the  node  voltage  to  the  branch  with  the  lowest 
i-mpedance  factor  and  therefore  tends  to  stabilize  the  node  voltage  against 
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numerical  oscillations.   Because  of  this  stabilizing  effect  this  last  criterion 

for  the  calc\ilation  of  the  node  voltages  seems  to  he  the  best.   However^  no 

detailed  analysis  comparing  the  last  two  methods  has  yet  been  made. 


3.   CALCULATION  OF  IMPEDANCE  FACTORS 

In  the  above  discussion  the  knowledge  of  the  values  of  Z  ,  Z  and  Z 

c"^   R      L 

for  each  branch  was  assumed.   Since  there  are  three  parameters  which  must  be 
defined,  we  need  three  sets  of  equations  for  each  branch.   This  makes  impossible 
the  redefining  of  a  set  of  impedances  during  each  step  of  the  analysis.   Though 
this  would  be  desirable,  it  is  not  absolutely  necessary.   The  initially  defined 
impedances  will  never  change  if  all  the  elements  are  linear  and  will  be 
relatively  constant  over  small  ranges  of  voltages  even  if  the  elements  are  non- 
linear.  In  the  latter  case  the  impedance  factors  must  be  redefined  when  the 
actual  branch  voltage  changes  start  to  differ  appreciably  from  the  ones  predicted 
with  this  set  of  impedance  factors.   If  we  adhere  to  the  restrictions  in  the 
topology  previously  mentioned  there  are  two  methods  which  can  be  used  to 
calculate  a  new  set  of  impedance  factors. 

Let  us  assume  that  we  have  the  condition 

or 

Under  this  assumption  it  may  be   recalled  that  we  will  have  for  each  branch 


A^Q..    =  A^Q..    =  A  0.. 
A  Ji  A  ji  A  ji 


Therefore 


Av"!".At. 

-7^-^^Ci*hi^W'^i  (3.1) 
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After  a  complete  numerical  cycle  all  the  terms  on  the  left  are  known  and 
therefore  the  impedance  factor  Zj  can  be  calculated. 

To  separate  the  three  impedance  factors  two  more  equations  are  needed. 
To  obtain  these  we  must  first  define  an  equivalent  incremental  capacitance, 
resistance  and  inductance  which  are  in  series  between  point  Qj  and  the  reference 
node : 


AQ.. 

^^  (3  =  2) 


^  CEji 


Av 

R^.  =-^  (3»3) 

Ji 

At. 
1 


1^=^^  (3A) 

(At.)2 


Combining  these  with  the  definition  of  Z  .,  Z^^.  and  Z       we  get 


At. 


^RJ  =  %  <3-'^ 


Z,  •  =  ^  '3.7) 


For  small  changes  in  voltages,  the  equivalent  elements  C^  ,  R^  and  L^  will  be 
constant ,      Hence 


Z^.  ^   At.  (3-8) 

Cj     1 
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Z   =  constant  (3»9) 


^Lj "  i:  (3.10) 


From  this  we  see  that  if  we  repeat  the  numerical  cycle  three  times  with  three 

different  values  of  time  increments  k-|At.,  k^At.  and  k^At .  we  will  be  able  to 

1   1   2   1      3  1 

get  three  values  for  each  Z..   By  using  the  variation  of  the  impedance  factors 

J 

rfith  At.    noted  in  equations   3-8  to   3«10^    we  are  able  to  write   three   equations 

Z 

k  Z^.    +  Z^  — -^  =  Z.  m  =  1,    2,    3  (3oll) 

m  Cj  R.    +     k  jm  y      >    -J  \->        J 


//here  Z.      is  the   impedance   factor  found  by  using  time   interval  k  At..      When 
jm  mi 


<:„  =  1.0,    the   solution  of  these  equations   yields 


k  k 
Z.    =  -  — -^ 


L       k^   -   k3 


"2    -   \        VA' 
.   1-k^-      k3-l 


(3^2) 


^c  =  -^^k^  (3a3) 


Z^=   \   -^L-^c  (3.1V 


The  condition  AQ  .  »  A.  Q,  .  can  be  realized  by  injecting  into  the 
.op  node  a  pulse  which  is  much  higher  than  the  normal  signal.   Since  this  will 
i'esult  in  large  changes  in  voltages,  all  nonlinear  elements  must  be  temporarily 
•eplaced  by  a  linear  equivalent  whose  parameter  values  are  the  same  as  the 
ncremental  parameters  of  the  nonlinear  element  during  the  last  normal  cycle » 
tecause  of  these  much  higher  than  normal  responses,  one  of  these  impedance 
actor  calculation  cycles  cannot  be  used  as  a  normal  cycle  in  which  time  is 
ncremented. 
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The  second  method  of  determining  the  impedance  factors  also  depends 
on  the  assumed  topology.   If  the  circuit  conforms  to  this  topology,  the  various 
impedance  factors  are  nothing  more  than  the  equivalent  impedances  looking  from 
the  point  in  question  toward  the  reference  node,,   If  we  recall  that  a  node 
impedance  factor  Z  was  already  defined,  we  can  use  this  formula  to  calculate 
three  impedances  for  the  th_ree  time  increments  k-j^At^,  -^t.    and  k  At^,   We  can 
then  again  use   equations  3.12  to  3„1^  to  calculate  an  equivalent  capacitance 
C  J    resistance  R^  and  inductance  L^  for  the  node„   Since  these  are  in  series, 
they  can  be  added  to  the  corresponding  elements  in  the  branch  above  this  node. 
In  this  way  all  the  impedance  factors  may  be  calculated  in  a  more  direct  manner 
than  that  of  the  first  method.   Since  this  last  method  is  also  much  easier  to 
program  and  much  more  conservative  of  computation  time,  in  what  follows  It  will 
be  assumed  that  this  latter  method  is  being  used  to  calculate  the  impedance 

factors. 

At  this  point  we  might  also  introduce  the  conditions  required  to 
change  the  time  increment  during  the  process.   We  must  impose  the  condition 
that  during  the  time  increment  change  the  voltages  across  the  three  impedance 
factor  elements  defined  in  equations  3  =  2  to  3,^4-  remain  constant; 

v,^  .  =  C,^  .  Q .  =  constant 


AQ. 

V  _  .  =  R^ .  —77-  =  constant 
PZ  j    ^  j  At 


2 

v^„.  =  Lv- .  — ~^  =  constant 
^J    ^J  (At)2 

Since  C   ,  R_  and  L^  are  also  constant  during  the  increment  size  change,  we 
conclude  that  the  terms 
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AQ.  A^Q. 

Q  ,      -^      and      ^ 

J        ^^  (At)^ 


must  also  remain  constant o   This  gives  us  our  desired  information, 


k.      VOLTAGE- CONTROLLED  SUBNODES 

In  the  previous  formulation  the  current- controlled  nonlinear  element 
has  a  decided  advantage  since  all  the  currents  are  determined  before  the  branch 
and  node  voltages  are  csuLculatedo   In  order  to  allow  voltage -controlled  devices 
to  be  handled  just  as  efficiently  it  is  desirable  to  formulate  a  similar  process 
which  will  be  compatible  with  the  above  but  in  which  the  role  of  the  current  and 
the  voltage  are  interchanged.   In  this  case  we  will  allow  KCL  to  be  violated. 

In  the  development  of  the  relevant  formulas  we  will  assume  that  we  are 
working  with  only  one  branch  extracted  from  the  circuit.   To  be  as  general  as 
possible,  we  can  postulate  that  this  branch  is  made  up  of  a  series  string  of  RLC 
elements  as  shown  in  Fig,  k,l.      The  current  generator  beside  each  parallel  unit 
represents  the  error  current  for  that  unit.   Since  the  nodes  between  the  parallel 
units  are  associated  only  with  branch  j  and  do  not  have  "levels"  assigned  to 
them,  we  will  call  them  'voltage -controlled  subnodes"  to  distinguish  them  from 
the  normal  nodes.   To  derive  the  necessary  relations  we  must  turn  to  the  dual 
concept  in  which  the  flux  is  defined  by 

cp  =  r  v(t)dt  (i^.l) 

0 

Then  for  linear  elements  we  have 


cp  =  LI 


^  =   RI  (i^.2) 

dt 


2 
dfcp  ^  1   -. 

dt^"^ 
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Figure  1^,1.   A  Branch  Consisting  of  M  Voltage-Controlled  Subnodes, 
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If  we  go  to  the  difference  notation  we  get  for  the  mth  parallel  unit 


L 
.  „    in 
cp  .  =  AQ  .  — — 
mi     mi  At. 
1 


ACp  .  =  AQ  .R  (4„3) 

mi     ml  m 


,2.     .„   ^*1 


^  V  =  '"Sai  -C" 


m 


where 


Acp  .  =  cp  .  -  cp  .    , 
mi    mi    m(,i-l; 


^2, 


A^cp  .  =  Acp  .  ~  Acp  ..  ^x  {h.h) 

mi     mi     m(i-l; 


and 


^22 
A^cp  .  =  A'^cp  .  -  A  cp 

mi      mi      m(,i-l; 


These  equations  suggest  that  we  define  a  set  of  "admittance  factors"  in  a 
similar  manner  as  the  impedance  factors  were  defined: 


PL  cp    . 

m  rai 


f         =  ^^  (4.6) 

PR         Acp   . 
m  mi 


AQ, 
Y. 


m       A  cp    . 
mi 


'Cmi  (i^„7) 
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^Pm  ~  ^PLm  ^  ^PRm  +  ^PCm  (^"^^ 


where  AQ^,   AQ^,  AQ^  are  the  charges  flowing  through  the  parallel  L^  R  and  C 
elements  respectively.   Also  define 


A^Q  .  =  AQ..  -  AQ  .  (i+  Q) 


where 


AQ   =  AQ   .  +  AO^  .  +  AQ^  . 
mi     Xjini    ^mi     Cmi 

2 
A  Q   corresponds  to  Av . .  in  the  main  method. 

Let  us  suppose  that  during  the  time  interval  (i-l)  an  error  of 

2 
^  "^mfi-l)  ^^sulted  at  subnode  m.   The  condition  for  the  compensation  of  this 

error  in  this  step  requires 


^2,  ,    , 

•1 


The  amount  of  flux  A^cp  .  required  to  do  this  is  given  by 


3       1 


Pm 


_^  ^m(i-l)  -  ^PLm(%(i-l))  -  (YpLm  ^  ^PRm^^^fi 


(i-l) 


(^ai) 


This  correction  formula  is  analogous  to  that  of  equation  21-6  and  hence  can  be 
treated  with  the  same  type  of  numerical  procedure.   Furthermore;,  any  future 
comments  on  equation  2.16  can  therefore  also  be  applied  with  minor  changes  to 
equation  ^,11.   Further  derivations  parallel  to  those  of  sections  2.2  and  2„3 
could  be  continued  here  which  would  then  make  it  possible  to  work  directly  with 
such  strings  of  parallel  elements  as  shown  in  Fig.  4.1.   However,  in  our  method 
we  can  be  content  with  just  being  able  to  handle  one  set  of  parallel  elements. 
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The  subscript,  m  will  be  retained  to  distinguish  charge  in  the  voltage -controlled 
subnode  from  the  other  charges. 

Besides  equation  i4-.ll  we  also  need  a  method  of  interconnecting  the 
voltage -controlled  node  with  the  rest  of  the  nodes  of  the  circuit.   To  obtain 
this  let  us  suppose  that  branch  j  which  contains  the  voltage- controlled  subnode 
has  a  set  of  impedance  factors  Z^.,    Z^.  and  Z^^  defined  for  it  and  is  regarded 
by  the  charge  distribution  method  of  section  2.2  as  a  normal  branch.   After  the 
charge  distribution  cycle  a  charge  AQ   has  been  assigned  to  flow  through 
branch  j .   The  increase  in  charge 


^  Si  =  ^Si  -  ^S(i-i) 


can  therefore  be  treated  exactly  the  same  as  the  error  charge  A  Q^^_-l)  and 
hence  the  two  can  be  summed  to  form  a  total  effective  error  charge  of 


^Sn(l 


P  2 

1-1)      ^m(x-l)      Ji 


(i^.l2) 


Hence  the  increase  in  <P  is  determined  by 


a3cp. 


1    .2 


mi   Y. 


Pm 


[4«»(1-1)  -  ^Pl(^^(1-1)'  -  (^PL  *   ^PR)^\(i-l), 


(l».13) 


Since 


Acp 


V 


mi 


mi    At, 


we  get  the  prediction  formula 


^Sii(i-l)   ^PLm  ^^ 
mi    m(i-l)      m(i-lj     Y.^At 


Pm  i 


^Pm  '"(l-^' 


Y    +  Y 
PLm    PRm 

^Pm 


Av 


m(i-l) 
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where 


%i  =  \±   -   \il-l)  {h.iy 


With  this  formula  the  voltage  across  the  parallel  elements  can  be  calculated 
and  consequently  the  currents  through  them  can  be  obtained.   The  total  charge 
%i  ^^"^  therefore  be  computed  and  the  error  charge  for  this  cycle  becomes 


A  Q^.  =  AQ   -  AQ  (^^,6) 


After  these  calculations  have  been  performed  for  all  voltage-controlled  subnodes, 

all  branch  voltages  are  known  and  therefore  the  node  voltage  calculation  part 

of  the  numerical  cycle  can  be  started. 

From  the  above  it  should  be  apparent  that  branch  j  which  contains  the 
voltage-controlled  subnode  can  also  have  the  full  complement  of  series  elements 

in  addition  to  this  subnode.   This  results  because  we  essentially  have  two 
error-correcting  mechanisms  associated  with  branch  j  and  therefore  it  can  have 
two  sets  of  elements,  one  in  series  and  another  in  parallel. 

By  the  use  of  this  voltage-controlled  subnode  technique,  restriction 
two  of  section  2.1  (only  one  branch  can  be  connected  between  two  nodes  except 
if  one  of  the  nodes  is  the  reference)  is  eliminated  since  we  now  have  a  method 
for  dealing  with  parallel  branches  without  violating  the  first  restriction.  Of 
course  this  technique  is  also  useful  for  voltage-controlled  elements  since  the 
voltage  across  the  elements  is  first  predicted  and  the  currents  are  calculated 
later. 


5,   DISCUSSION  OF  CORRECTION  FORMULA 

In  section  2.2  we  derived  equation  2,l6  which  will  compensate  for  a 
residual  error  between  the  branch  and  node  voltages  from  the  previous  step. 
This  formula  was  developed  so  that  a  circuit  of  the  assiomed  topology  and  with 
correctly  defined  impedance  factors  will  respond  so  as  to  have  a  zero  resiilting 
error  between  v  .  and  vt. .   If  the  impedance  factors  are  not  completely  correct 
due  to  the  presence  of  nonlinearities,  equation  2.l6  will  not  make  the  proper 
correction. 

5„1  Error  Term  for  a  Branch 

To  find  the  factors  related  to  this  correction  formula,  let  us  examine 
the  error  that  will  be  produced  when  the  impedance  factors  are  not  correctly 
defined.   To  be  able  to  do  this  we  make  the  following  assumptions: 

(1)  The  branch  considered  has  impedance  factors  such  that 
V  is  relatively  independent  of  v  . 

(2)  Z',  the  impedance  factor  actually  used,  and  Z,  the 
correct  impedance  factor,  are  related  by 

Z'  =  Z  +  AZ 

where  AZ  is  the  error  in  Z' .   Similar  definitions  hold 
for 

^'C    \'   ^C  ^R'  ^R'  ^R'  ^'  ^L  ^^  ^ 

2 
(X)      Q  ,,  AQ.  -,  and  A  Q.  ,  are  known. 

^~"  1-1-^      1-1  1--L 

iS)     The  error  voltage  after  the  last  step  in  the  numerical 
process  was 
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(5)   The  node  voltage  v  will  be  increased  by  an  amount 

Av,  .. 
ki 

The  IQ  method  will  then  produce  the  terms 


and 


4«i  -  f^  hl-l)  -  ^C  ^  -  (^H  ^%^^]  (5.1) 


Av  At 


a3q'  =  a3q^  +  a3q^  (5„3) 


The  correct  third  difference  of  charge  should  be 


^%=fK-X-c%^-(^C-K)%i|  (5..) 


4^  '-   -^  (5.5) 


The  error  is 


a3q.  =  a3q.  -  a\'  (5^6) 


By  the  substitution  of  appropriate  definitions  this  b 


ecomes 


3        AZ  ^^n  -  ^r^  Z^  -  Z^T   o 

^\  =  Z(^Taz)  (^\-1  -   ^\i)^t  .  -^xfr^  ^^i-l  -  Z(Z.AZ)   ^^^i.-l 

(5.7) 
The  new  error  voltage  for  step  (l)  is  given  by 
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Av.  =  V,  .  -  V. 
X    kx    1 


By  using  previously  derived  relations  we  get 


Av^  =  Z^  (5.8) 


Hence 


ZAZ^  -  Z  AZ  AQ.  ,    Z^AZ  -  ZAZ^  A  Q.  -,         ' 
AZ    /,       ,    N      C    (J     x-1  ,  _L L   ^x-1 

^^i  =  zVaZ  ^^^-1  ^  ^""ki^  ^   Z  -kAZ      aT  "■    Z  +  AZ     At  ; 

(5.9)   ;| 

The  first  term  indicates  that  if  AZ  is  small,  the  error  due  to  i^^^^j^  +  ^^ki^ 
will  be  negligible.   The  second  and  third  terms  are  directly  related  to  the 
current  and  the  current  increase  respectively  in  the  branch.   It  may  be  obserA-ed  '1 
that  these  can  have  much  larger  magnitudes  than  the  first  term.   Nevertheless 
they  do  exhibit  the  convenient  property  of  being  slowly  varying  functions  of 
time.   This  is  a  consequence  of  assuming  that  the  impedance  factors  are  rela- 
tively constant  even  thovigh  a  nonlinear  element  may  be  present  and  that  the 
current  will  normally  not  vary  greatly  from  one  numerical  step  to  the  nexto 

5.2  Compensation  of  the  Error 

The  error  caused  by  the  last  two  terms  in  equation  5-9  is  mostly 
attributable  to  the  fact  that  the  impedance  factors  were  developed  so  that  the 
circuit  responds  accurately  to  a  pulse  of  width  At  which  for  linear  elements 
connected  in  the  assumed  topology  may  be  of  any  height.   If  a  step  function  is 
applied,  the  circuit  will  respond  accurately  initially,  but  as  the  effects  of 
the  transient  decrease,  new  impedance  factors  corresponding  to  longer  At's  should 
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be  used  and  eventually  when  the  circuit  is  again  in  a  steady  state  the 
impedance  factors  should  merely  he  defined  from  the  resistances  of  the  circuit. 

We  note  that  as  the  effects  of  the  transient  are  dying  out  the  cur- 
rents will  become  more  and  more  constant  from  one  step  to  another „  Therefore 
the  terms 

ZAZ^  -  Z^  AQ._^   Z^AZ  -  ZAZ^  A^Q. 
Z  +  AZ     aT  +    Z  +  AZ     A^ 

may  be  compensated  by  the  introduction  of  a  slowly  time-varying  factor  into  the 
correction  formula.   Since  the  error  impedance  factors  of  the  form  AZ  are  not 
known  by  the  numerical  process,  this  factor  would  have  to  be  adjusted  by 
Observing  the  trend  of  the  error  voltage  Av.  and  making  changes  to  minimize  Av 
in  the  next  few  steps.   The  introduction  of  this  adaptive  variable  eliminates 
the  need  of  a  multiplicity  of  impedance  factors  which  would  be  impractical  to 
use. 

In  the  correction  formula  2.l6  we  replace  the  term  Av,    ,  bv 

(i-1)  -^ 

^""(i-l)  +  \   ''^^^^  \   is  the  new  factor  mentioned  above.   Whenever  the  circuit 
is  in  a  steady  state,  then  at  each  branch 


A^Q.  »  0 

1 


Av.  ~  0 

1 


This  requires  that 


""a  "  ^C  ^  (5.10) 


When  the  transient  analysis  is  started  this  condition  must  be  forced  and  during 
transients  v^  should  be  adapted  so  that  when  a  steady  state  is  reached  this 
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condition  will  again  hold.   By  the  introduction  of  this  correction  factor  the 
response  of  the  circuit  to  the  transients  is  in  no  way  impaired;  the  impedance 
factors  defined  for  At  look  after  the  fast  transients  and  v^  adjusts  the  circuit 
to  the  slower  transients. 

It  should  be  emphasized  at  this  point  that  the  error  voltage  Av^  at 
each  branch  is  known  after  each  numerical  step  and  hence  can  be  controlled „   If 
an  error  at  a  branch  is  too  large,  one  can  always  repeat  the  step  using  a 

smaller  time  increment.   Since  Z  varies  directly  with  At  and  since  the  current 

AQ. 

— -   remains  constant  during  the  step  size  change  the  error  term 

At 

AQ. 
^C  At 

will  also  decrease  with  a  decrease  in  the  time  interval.   Hence  the  introduction 

of  V  is  not  absolutely  necessary,  but  its  use  will  increase  the  speed  of 
a 

computation  considerably. 


5,3  A  Predictor-Corrector  Type  Method 

To  improve  further  the  properties  of  the  IQ  method  we  can  introduce 

a  predictor -corrector  type  procedure  which  is  somewhat  analogous  to  those  used 

[6] 


in  the  integration  of  differential  equations 
eqixation  2„l6  as 


We  redefine  the  correction 


4^--f 


Av 


m 
Ei 


'C  At 


^\  ^  ^r) 


At 


m  =  1,  2,  00,,  M 


(5.11) 


where  M  is  the  number  of  times  the  process  described  in  section  2  is  repeated 
for  each  increment  in  time.   The  procedure  is  to  calculate  a  value  for  Av^^ 
for  each  node,  use  this  in  the  distribution  of  the  currents  and  calculate  all 
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the  node  voltages.   An  error  voltage  Av.  will  result  which  is  then  used  to  define 

2 
a  second  value  ^v^^.   This  cycle  is  then  repeated  M  times^,  the  last  "being  the 

one  which  will  he  accepted  as  yielding  the  charges  and  voltages  for  this  time 
increment . 

The  problem  is  to  determine  a  set  of  expressions  for  the  terms  Av"^  . 

Ei 

This,  in  the  general  case,  could  become  quite  complicated  but  the  results  could 

prove  to  be  valuable  here  in  the  same  way  as  those  obtained  by  including  more 

terms  in  the  predictor  formulas  in  the  integration  of  differential  equations. 

To  simplify  matters  we  will  restrict  the  discussion  to  the  case  in  which  M  =  2, 

a  value  which  would  probably  be  close  to  the  optimum  when  computation  time 

versus  decrease  in  error  is  considered.   The  expression  for  Av„.  then  has  the 

El 

form 


Av^.  =  V  +  a  Av.  -, 
El    a    1  i-1 


a^  can  be  either  a  constant  for  the  branch  or  an  adaptable  factor  which  can  be 
slowly  adjusted  during  the  analysis  so  that 


Av!  ~  0 
1 


2 
The  expression  for  Av„.  becomes 

El 

Av2.  =  ^^  ^  a^Av._^  +  a^Av^ 

where  a^  is  again  a  constant  or  adaptable  factor. 

Since  the  three  variable  v  ,  a  and  a  all  are  directed  towards 
decreasing  the  error  Av^,  it  is  obvious  that  they  are  highly  interrelated. 
But,  the  role  that  they  play  is  different:   v  is  related  to  the  algebraic  sum 
of  the  error  over  a  number  of  time  intervals  and  a  and  a  are  related  to  the 
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numerical  oscillations  and  damping.   Since  these  parameters  must  be  varied 
slowly,  it  is  apparent  that  the  numerical  process  must  be  slow  enough  to  allow 
the  adaptation  of  these  parameters  to  the  response. 

This  predictor-corrector  type  procedure  was  not  studied  further,  but 
an  analysis  of  the  factors  introduced  and  their  relation  to  stability  would 
probably  prove  useful. 


6c   ERRORS  AND  STABILITY 

6.1  Circuit  Error  and  Response  Error 

Up  to  the  present  we  have  been  talking  about  an  error  A^oltage  Av. 
at  each  branch  or  an  error  current  A  Q./At  at  a  voltage-controlled  subnode. 
For  purposes  of  discussion  we  will  call  these  the  "circuit"  errors »   There  is 
yet  another  type  of  error  which  we  will  name  the  "response"  error,,   The  latter 
refers  to  the  difference  between  the  theoretical  value  of  a  current  or  a  voltage 
at  a  given  time  and  the  value  of  the  same  current  or  voltage  obtained  by  the 
numerical  process  at  the  same  time  instance. 

Previously  it  was  shown  that  the  circuit  error  can  be  controlled  by 

specifying  a  maximum  allowable  error  between  any  branch  and  the  node  to  which 

it  is  connected.   If  we  have  K  nodes  connected  in  the  worst  case  obtainable^ 

namely  all  in  series,  and  if  at  each  node  there  is  an  error  Av    and  since  the 

max 

voltage  at  the  branch  connected  to  the  highest  level  node  must  be  within  +  Av   , 

—   max 

we  can  conclude  that  the  circuit  error  of  any  node  voltage  is  less  than  —  Av 

2   max 

Normally  this  error  would  be  much  smaller  and  in  the  specific  case  the  exact 

upper  bound  can  be  determined  from  the  topology  of  the  circuit.   We  can  also 

conclude  that  whenever  the  circuit  is  in  a  steady  state  the  response  error  will 

be  less  than  this  maximum  circuit  error. 

During  transients  the  last  statement  does  not  hold,  but  we  can  define 

a  rough  estimate  of  the  upper  bound  for  the  response  error  during  the  transient 

by  considering  the  special  case  in  which  a  linear  capacitance  is  one  of  the 

branches  below  the  node  whose  response  error  we  wish  to  find.   Suppose  that 

the  theoretical  current  flowing  through  the  capacitance  is  given  by  l(t)  and 

AQ^ 
that  the  current  calculated  by  the  IQ  method  is  — r-  where 

iAt  <  t  <  (i  +  l)At 

Also  define 
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-31^- 


1 


Then  the  error  in  the  charge  on  the  capacitance  is 


^^  m=l      m=l 


If  we  assume  Av.  =  0  (ioe,,  v   =  v.),  the  response  error  becomes 

1  i^l        1 


1  .  „    1  ^  .2, 


^^E1  =  C^\  =  C  ^AS 


-,  ^  m 

m=l 


But  the  maximum  circuit  error  possible  is  Av    and  consequently 

max 


2 
^^i 
C       max 


If  there  is  no  cancellation  of  the  circuit  error  from  one  time  interval  to 
another,  we  get  the  maximum  response  error 


Av^.    =  iAv 
Eimax      max 


In  practical  situations  this  would  be  much  lower  and  fiirthermore  a  much  closer 

upper  bound  for  a  specific  circuit  can  be  calculated  by  observing  the  circuit 

error  at  each  time  interval  during  the  transient.   From  this  it  is  also  apparent 

that  the  adaptive  factors  v  and  a^  should  be  adjusted  so  that  small  numerical 
^  a     1 

oscillations  will  result.   This  will  cause  the  circuit  error  to  oscillate 
between  positive  and  negative  values  and  therefore  produce  a  much  smaller 
response  error.   It  is  also  felt  that  a  much  lower  upper  bound  for  the  response 
error  could  be  arrived  at  if  the  energy  conservation  criterion  is  used  to  cal- 
culate the  node  voltage,  but  this  was  not  examined. 
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6,2  stability 

The  stability  of  the  IQ  method  has  not  yet  been  examined  in  detail, 
but  from  the  discussion  of  errors  we  can  make  a  few  observations  about 
stability.   First  let  us  recall  that  in  the  normal  predictor-corrector  methods 
of  numerical  analysis  there  is  available  an  error  term  which  is  used  to  control 
the  increment.   Consider  the  example 

Predictor     x   .    =  x    ,      ^^    +   Atx  ,    ^ 
Pi    c(i-l)      c(i-l) 

Corrector     x    .    =  x    ,      , .    +  —   (x       +  x    ,  ) 

ci    c(i-l)    2  ^V  ^  ''p(i-l)^ 

where  the  term  (x^.  -  x^.)  is  the  error.   It  is  to  be  noted  that  this  is  a 
relative  error  since  the  corrected  as  well  as  the  predicted  value  may  be 
erroneous.   When  both  terms  are  in  error  but  the  difference  between  them  is 
small,  the  solution  will  possibly  become  numerically  unstable. 

On  the  other  hand,  in  the  IQ  method  the  error  which  controls  the  step 
size  is  the  circuit  error,  an  absolute  error.   We  have  already  seen  that  we  can 
set  an  upper  bound  for  this  error.   Since  we  can  assume  the  circuit  error  is 
initially  within  the  allowed  bounds,  then  for  each  succeeding  step  the  error 
will  also  be  within  these  bounds.   The  last  statement  is  justified  if  we  assume 
that  there  are  no  discontinuities  in  any  v-i  characteristics  and  if  we  recall 
that  a  decrease  in  the  time  increment  implies  a  decrease  in  the  circuit  error. 
Therefore  as  long  as  we  have  continuous  functions  for  the  v-i  characteristics 
the  solution  will  not  be  able  to  "blow  up." 

The  problem  of  instabilities  in  the  integration  of  the  differential 
equations  of  a  circuit  are  usually  related  to  the  minimum  time  constants  of 
the  circuits  which  will  cause  instabilities  if  the  time  increment  is  too 
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largeo^"'"-'  For  the  IQ  method  this  does  not  appear  to  he  a  problem.   Consider 
for  example  the  tranch  in  Fig.  6.1  which  has  a  very  low  series  inductance  and 
a  high  resistance: 


L. 
Lj   At 


Rj    J 


Z  ~  0 


Figure  6.1. 

If  the  node  voltage  varies  slowly  with  time  and  a  large  At  is  heing  used  then 
Z  is  insignificant  compared  to  Z^^  and  consequently  Z^  is  all  but  ignored  by 
the  method.   If  a  large  transient  (with  a  finite  slope)  appears  at  the  node. 
At  will  be  decreased  until  the  node  voltage  change  is  only  incremental.   Now  Z^ 
has  a  much  higher  value  and  might  even  dominate  Z^.   As  the  transient  disappears. 
At  is  again  increased  in  size  and  Z^   again  assumes  an  insignificant  role.   Thus 
we  see  that  the  small  inductance  will  only  have  an  effect  on  the  circuit  during 
fast  transients. 


7.   RELAXATION  OF  TOPOLOGICAL  RESTRICTIONS 

Suppose  that  it  is  desired  to  perform  an  analysis  of  a  bridge  type 
circuit  shown  in  Fig.  7»la.   The  elements  marked  L  are  assumed  to  be  lower 
impedances  (e.g.,  10  ohms)  than  those  marked  H  (1000  ohms)o   Two  possible 
redrawings  of  the  circuit  according  to  the  method  described  appear  in  Figs.  7„lb 
and  7»lc.   Both  forms  of  the  circuit  can  be  analyzed  by  the  IQ  method  but  for 
fast  transients  considerably  more  difficulty  will  be  experience  with  Fig„  7olc 
than  with  Fig.  7.1b.   It  may  be  recalled  that  the  impedance  factors  were 
derived  from  the  "down-directed"  equivalent  impedances.   The  impedance  factor 
at  the  top  node  would  therefore  be  approximately  kO   ohms  for  Fig.  7 .lb  and 
about  500  ohms  for  Fig.  7.1c,  whereas  the  correct  impedance  is  a  little  less 
than  kO   ohms.   If  a  voltage  step  is  applied  to  the  circuit  of  Fig.  7olc,  the 
injected  current  will  be  grossly  underestimated  and  the  process  of  establishing 
equilibrium  would  take  much  longer  than  in  the  case  of  Fig.  7.1b  where  it  would 
be  established  in  very  few  numerical  cycles. 

The  conclusion  from  the  above  is  that  if  a  circuit  does  not  have  the 
desired  topology  a  wise  rearrangement  of  the  circuit  will  possibly  make  it 
amenable  to  as  fast  a  solution  as  that  of  a  circuit  which  has  the  restricted 
topology.   There  seem  to  be  a  few  general  rules  which  will  usually  lead  to 
convenient  configurations  if  the  original  circuit  does  not  have  the  desired 
topology: 

(1)  If  possible  the  highest  level  node  influencing  the 
portion  of  the  circuit  with  the  undesirable  topology 
should  have  as  low  an  impedance  path  to  the  reference 
as  obtainable . 

(2)  If  possible  the  circuit  should  be  arranged  so  that  the 
lowest  impedance  branches  lead  to  the  lowest  possible 
node  levels . 
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a)  ORIGINAL 


k3 

L  <  R. 


L<R. 


R3<H 


L<^K 


Hi 


L>R, 


b)  CORRECT  REDRAWING 


lSR, 


L>. 


R^^H 


C)  WRONG  REDRAWING 


Figure  7„lo      Example   of  a  Circuit  Violating  the  Assumed  Topology. 
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By  allowing  these  variations  in  the  topology,  we  must  also  restrict  the 
forcing  function  voltage  to  incremental  changes  from  one  time  interval  to  the 
next  since  the  impedance  factors  will  no  longer  give  the  exact  response  to  a 
high  pulse.   In  a  specific  circuit  the  amount  of  difficulty  the  method  will 
encounter  can  be  predicted  since  it  is  directly  related  to  how  well  the 
impedance  factors  approximate  the  actual  equivalent  impedances. 


8o      NONLINEAR  AND  TIME- VARYING  PARAJvETERS 

8.1  Soiirces 

A  voltage  source  is  already  Included  in  the  IQ  method  and  a  constant 
current  source  does  not  present  any  difficiilties  since  this  requires  only  that 
a  charge  I^At  be  subtracted  from  one  node  and  added  to  another  before  the 
start  of  the  charge  distribution  cycle.  With  constant  sources  there  are  no 
problems  (except  possibly  during  the  initialization)  but  if  they  are  time- 
dependent  or  proportional  to  currents  or  voltages  in  another  part  of  the  circuit 
the  method  might  have  trouble.   If  a  voltage  soaarce,  for  instance,  is  time- 
dependent  and  during  one  numerical  cycle  experiences  a  large  change,  this 
change  will  not  be  taken  into  accoiint  by  the  first  charge  distribution  cycle 
since  no  impedance  factors  are  changed.  When  the  node  voltages  are  then  cal- 
ciilated  there  will  be  a  large  error  at  the  branch  with  the  voltage  source. 
Since  this  error  will  then  be  diminished  by  the  next  numerical  cycle  through 
the  error-correction  mechanism,  we  see  that  the  response  of  the  circuit  to  all 
varying  voltages  will  be  delayed  one  cycle,   A  current  source  will  also  produce 
the  same  result,  but  with  the  slight  advemtage  that  the  node  to  which  it  is 
connected  will  respond  immediately  and  the  delay  will  appear  in  a  higher 
level  node.   This  advantage  can  often  be  put  to  good  use.   The  rule  in  general 
seems  to  be  that  only  incremental  changes  should  be  allowed  from  one  time 
interval  to  another  for  varying  sources, 

8.2  Nonlinear  Elements 


In  general  we  can  specify  the  characteristics  of  a  nonlinear  element 


by 


-^^  =  f (v)    or    V  =  g  \-^^  m  =  0,  1  or  2        (8,l) 

(Atr  \(Atr/ 


-ko- 
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The  second  form  is  the  one  which  falls  naturally  into  the  domain  of  the  IQ 
method  but  unfortunately  it  is  not  always  the  most  convenient  to  use„  For 
instance,  in  a  tunnel  diode  the  relation 

is  multivalued.   In  such  cases  the  voltage-controlled  node  formulas  developed 
earlier  can  be  used  in  conjunction  with  the  first  of  the  above  form,ulas  which 
is  in  this  case  single  valued.   The  important  thing  is  that  the  prediction  and 
circuit-error  correction  mechanism  is  supplied  by  the  method,   therefore 
avoiding  the  need  of  a  repetitious  process  such  as  Newton's  method  to  determine 
the  independent  variable  when  the  dependent  variable  is  given.   This  can  prove 
to  be  an  important  time-saving  factor  during  computation. 

80 3  Semiconductor  Diodes 

Since  the  diode  and  transistor  are  the  most  important  nonlinear 
elements,  they  will  be  discussed  here  in  more  detail.   The  dc  characteristics 
of  an  ideal  diode  are  given  by 


l  =  lQ[e^-l]  A.  =  ^  (8.2) 


q 


The  Incremental  conductance  is  given  by 


V 


(8,3) 


^D   dv    X 
During  the  transients  the  incremental  diode  capacitance  given  by 

^  "  dv 


(8.i|) 
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also  becomes  important «  When  in  the  nonconducting  region  the  diode  exhibits 
a  depletion  layer  capacitance  of  the  form 


C  = -^ (8„5) 


where  v  is  the  built-in  potential,  n  is  the  grading  constant  which  is  usually 
D 

1    1 
equal  to  —  or  -  and  C-.  is  a  constant. 

The  diffusion  type  capacitance  when  in  the  conduction  region  is 

C  =^^==g  ^  (8.6) 

D   dv   ^D  dl  ^    ^ 

The  term  — ^  is  a  time  constant  which  can  usually  be  assumed  constant  as  a 
dl 

first  approximation  and  equal  to  the  diode  constant  T    Hence 

C  ~  g^T^  (8.7) 

Normally  only  one  of  the  two  capacitances  will  be  important  in  a  particular 
application.   Tn  switching  circuits  C  is  most  important  and  usually  C^.,  is  small 
enough  that  it  can  be  approximated  by  a  constant. 

If  T   cannot  be  considered  a  constant  the  relation  Q  =  f (v)  must  be 
obtained  by  niimerical  integration  during  the  transient: 

Cj^(v)dv=j   gj^(v)Tj^(v)dv  (8,8) 

0  0 


If  T   can  be  assiimed  constant  we  get  the  relation 


v 
Q  =  .,l/  (8.9) 
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When  the  depletion  layer  capacitance  parameters  v^  n  and  C  can  be  assumed 
constant,  and  if  the  diode  is  reverse  "biased,  we  get  the  formula 


Q  = 


n  -  n  *- 


(8a0) 


From  a  programming  point  of  view  the  expressions  for  1,    g  ,  C  and  Q  are  in  a 
very  convenient  form  since  they  are  so  highly  interrelated »   This  gives  the 
very  important  result  that  no  further  simplifying  assumptions  which  are  usually 
made  in  a  normal  analysis  will  be  required. 

To  complete  the  model  for  the  diode  we  can  include  the  leakage 
resistance  R^  and  if  the  reverse  breakdown  is  important,  a  current  of  the  form 


^B  = 


'0 


1  -  (^)- 


(8.11) 


where 


1 


1  -  (^)- 
^B 


is  the  multiplication  factor  and  v_  is  the  reverse  breakdown  voltage o 

Since  all  the  above  elements  are  in  parallel  and  since  the  current 
equation  is  in  the  form 


we  must  use  the  voltage-controlled  subnode  methods   The  relevant  parameters 
are  given  by 
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♦ *  9 


>  ^dS 


^D ^C,  V 


a)  IN  TERMS  OF  DIODE  PARAMETERS 


.      AO: 

6\     J 


b)  IN  TERMS  OF  ADMITTANCE  FACTORS 


Figure  8.1<,      Semiconductor  Diode  Model. 


c 

(8.12) 

y    —     ^  "V   — 

PC   At         P  ~  Sj-  +  ^  +  ^ 

c 

The  prediction  formula  h.lk   therefore  reduces  to 

4^(i-l)   ^PR 
-i  =  -i-1  -  ^-(i-1)  -  -Y^  -   -f  AV(  ._^)  (8„13) 

The  spreading  resistance  and  the  wiring  inductance  have  not  yet  been  included^ 
hut  these  present  no  prohlem  since  we  can  still  specify  a  complete  set  of 
series  elements  for  the  branch  containing  the  diode. 

Q.h     Transistors 

Of  the  different  large  signal  models  available  for  the  transistor j,  it 
was  found  that  for  the  IQ  method  the  Ebers  and  Moll  model  gives  the  best 
balance  between  accuracy  and  efficiency  in  computation.   For  this  reason  the 
technique  available  to  handle  this  model  will  be  developed  in  detail. 

The  equations  of  a  PNP  transistor  for  the  dc  case  are 

I       ^  T       ^ 

^  =  -r^(e-  -i)..^^^(e-  .1)       (8.ia) 

I        ^  T        ^ 

T      -  ry  EO    ,  kT     ^  s        CO    ,  kT     ,  .  /o  .  n 

^c  -  ^n  1  -a^^   (^    -  1)  -  r^^^  ^^    -  ^^       (Q"^5) 

The  notation  and  the  sign  convention  is  that  used  by  Ebers  and  Moll^^  with  the 
exception  that  the  symbol  9  is  replaced  by  v.   Therefore  these  will  not  be 
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discussed  here  in  detail.     For  future    convenience  let  us  make  the  following 
definitions  °, 


X  = 


kT 


-1 


■^ 


EO 


0      1  -  a 


^1 


Therefore 


E 


^=^o(-     -^' 


-I 


CO 


'CFO      1  -  ajXt-j- 


^CF  =   ^CFO^^        -   ^^ 


^  =  ^r  -  "-^' 


HF  I   CF 


-■-C        '^N-'-EF  "^   """CF 


(8a6) 


A  conductance  for  the  emitter  junction  and  one  for  the   collector  jiinction 
are  defined  as  follows  °, 


DTfl 


S 


E        1 5v, 


"SiF  "^  "S:: 


FO 


V  =const 
c 


(8,17) 


81, 


gn    = 


C       \5v. 


■'"CF  '^   "'"CFO 


Vw=  const 
ill 
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When  the  junctions  are  reverse  biased  the  leakage  resistance  R  and  R  for 

the  emitter  and  collector  junctions  respectively  become  important  and  therefore 

must  be  included  in  the  equivalent  circuit. 

During  transients  either  a  diffusion  capacitance  or  a  depletion  layer 
capacitance  come  into  play  at  each  junction,  depending  on  whether  the  junction 
is  forward  or  reverse  biased.   The  two  diffusion  capacitances  which  are  char- 
acterized by  two  time  constants  T^  and  T^  are  given  bv 

E      C  "^ 


C    =  T  a; 
DE    E^E 


^CE  ~  ^C^C 


(8.18) 


The  time  constant  T^  and  T^  can  be  assumed  constant  as  a  first  approximation 
and  therefore  similar  charge  versus  voltage  relations  as  equation  8„9  for 
diodes  can  be  obtained.   The  depletion  layer  capacitances  are  given  by 


o  -    °o= 


TE    /    v„  \n^ 


11-^ 


V    / 

DE/ 


(8,19) 
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where  C^^  and  C^^  are  constants  of  the  transistor,  v^^  and  v   are  the  built-in 
potentials  for  the  emitter  and  collector  junctions  and  n^  and  n  are  grading 
constants. 

Combining  all  of  the  aboA^e  we  have  the  model  shown  in  Fig,  8,2a , 
The  reason  for  the  unusual  orientation  will  become  apparent  later.   We  note 
two  important  properties  of  this  equivalent  circuit:   the  collector  and  the 
emitter  junctions  have  the  same  configurations  and  similar  elements  as  the 
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Figure   8.2.      IQ  Model   of  a  Transistor, 
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equivalent  circuit  for  the  diode  except  for  the  current  generator,  and  the 
circuit  is  in  the  exact  form  required  by  the  voltage-controlled  subnode  method 
of  section  k.      The  first  of  the  above  indicates  that  we  could  easily  use  the 
method  developed  for  diodes  in  section  8.3„   But,  the  second  property  suggests 
that  a  m.ore  general  technique  designed  specifically  for  transistors  can  be 
formulated „ 

Since  we  have  no  inductive  element  in  the  circuit,  equation  4,11 
reduces  to 


a\ 


mi   Y^ 
Pm 


Y 
PRm  „2. 


.4Sr.(l-l)J  -  tf  "^Vi-D  <^-'°' 


Pm 


To  put  this  in  terms  of  voltages  we  invoke  the  relations 


V 


At  ^^  =  -A^  ^^  =  ~^         (8-21) 


Therefore 

2 

^  Vi  -  -Y-aT-  -  ^7  ^%(i-l)  (8»22) 

Pm         P     ^    ^ 

For  the  transistor  we  have  two  parallel  units,  unit  1  at  the  collector  junction 
and  unit  2  at  the  emitter  junctiono  The  admittance  factors  for  these  two  units 
are  defined  by 

PRl   ^C   R 
c 

„      Sc   ^DC  /ON 

^PCl  =  -A^  -^  ^  (8-23) 


Y   =  Y    +  Y 
PI    PRl    PCI 
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Y    =  a:   +  — 
PR2    ^E    Rg 


PC2    At    At 


Y   =  Y    +  Y 
P2    PR2    PC2 


The  voltages  are  given  by 


a2 
V-,  .  =  V  .  =  V,  /  .  ,  N  +  Av-,  /  .  .X  +  A  V 
li    ex    l(i-l)     l(i-l)      li 


(8.25) 


^2i  ^  -^Ei  ^  ^2(1-1)  "^^2(i-l)  ^^^2i 


V  .  =  -v^ .  =  V^  /  .  .  \  +  Av„  /  .  ,  ^  +  A  V^_. 


The  two  current  generators  a   1        and  «  I   can  be  absorbed  into  the 
"error"  current  generators  A^Q^/At  and  A^Q^/At^   It  may  be  recalled  that  A  Q^ 
was  the  error  charge  as  defined  in  equation  4„9  and  AQ  was  the  sum  of  the 
charges  flowing  through  the  parallel  R,  L  and  C  elements.   For  the  transistor 
the  relations  become 


a\  =  AQ^^  -  AQ^ 


A^Q^  =  AQ^^  +  AQ^  -  AQ^ 


where  AQ   and  AO^^  are  defined  in  Fig,  8,2b  and  AQ^  and  AQ^  are  given  by 


v^  Av^       Av^ 


V  Av        Av 
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There  is  still  one  remaining  inconA/enience  which  concerns  timing., 
If  the  normal  IQ  procedure  is  followed,  the  current  will  first  be  distributed 
and  therefore  at  the  transistor  we  get  the  terms 


'^CTl  =  '^CT(l-l)  +  ^^%Tx 


V-^(l-l)  ^^\l 


The  error  term  which  is  to  be  used  in  equation  8„22  then  becomes 


^^li  =  ^'^Ki-D  -  ^'^cTi 


4^i=^%(i-l)  --^'Si-^^^Tl 


8.28) 


If  the  new  voltages  v^^  and  v^^  are  calculated  from  this,  there  is  a  possibility 
of  a  larger  error  in  AQ^  and  AQ^  since  the  correction  formula  does  not  take 
into  account  the  fact  that  OC^I^^t   and  aj^I^^^At  will  have  different  values  in 
time  interval  i  from  those  in  time  interval  (i-l)„   Consequently  the  two  cur- 
rent generators  will  always  be  lagging  behind  the  voltages  by  one  numerical 

cycle  during  transients „   To  compensate  for  this,  values  for  v,   and  v   can 

li      2i 

be  predicted  from  errors  in  Interval  (i-l)  and  these  may  then  be  used  to  predict 
A^alues  for  0!^I^.^t  and  Qj^IgpAt  for  the  time  interval  i„   The  changes  from 
interval  (i-l)  to  i  can  then  be  included  in  the  error  formulas  8o28o   It  is 
also  apparent  that  a  predictor-corrector  type  method  such  as  described  in 
section  5,3  could  be  used  here  to  further  reduce  this  numerical  lag  effects 

To  see  how  the  transistor  equivalent  is  incorporated  into  the  rest  of 
the  circuit  we  must  recall  the  assuro,ption  that  only  one  branch  of  the  main  cir- 
cuit is  involved  when  dealing  with  a  set  of  voltage-controlled  subnodes.  Hence 
the  entire  equivalent  for  the  transistor  except  the  base  is  included  in  one 


-52- 


branch.   The  collector  bulk  resistance   R   and  any  possible  lead  inductance 

cc 

can,  as  with  the  diode,  be  specified  in  the  series  elements  associated  with 
this  branch.   A  second  branch  is  required  for  the  base  and  therefore  it  would 
contain  the  base  resistance  R,-,  o   Since  the  voltage-controlled  subnodes  such 
as  that  in  the  center  of  the  transistor  model  are  normally  not  available  to 
the  main  IQ  method,  a  special  set  of  nodes,  one  per  transistor,  must  be 
defined  so  that  the  channeling  of  the  charge  which  flows  through  the  base 
branch  to  the  internal  transistor  subnode  can  be  accomplished. 

The  orientation  of  the  transistor  equivalent  shown  in  Fig„  8.2  is 
mainly  for  purposes  of  being  able  to  determine  the  approximately  correct 
impedance  factors  and  to  help  prevent  the  current  generator  lags  mentioned 
above.   The  parallel  admittances  of  the  transistor  equivalent  can  be  converted 
to  series  impedance  factors  by  the  method  described  in  section  3«   These  are 
illustrated  in  Fig.  8.3-   Normally  Z  would  be  fairly  small  and  capacitive 
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Figure  8.3 
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because  of  a  large  C^^,  and  Z^^  would  be  large  except  when  the  transistor  is 
in  saturation „   Thus  if  Z^   is  small  the  values  determined  by  the  impedance 
factor  calculation  method  of  section  3  will  be 

^2'\*h^   "CT  "  ^cT  (8.29) 

S'^l'\'  \o  (8.30) 

Because  of  the  normally  high  Aralue  of  Z^^  as  compared  with  R  ,  the  defined 
impedance  factors  are  a  good  approximation  of  the  actual  AQ  versus  Av  behavior 
at  the  transistor  terminals.   Therefore  this  orientation  of  the  transistor 
equivalent  will  normally  be  used„   As  an  example  of  how  a  transistor  circuit 
would  be  redrawn  before  analysis  by  the  IQ  method,  an  asymmetric  flipflop 
circuit  is  shown  in  Fig..  8„4, 

8„5  Transmission  Lines 

When  working  with  high-speed  circuitry  one  invariably  meets  a  situation 
in  which  a  signal  delay  is  involved.   Since  an  equivalent  transmission  line  will 
usually  represent  the  phenomenon  adequately,  a  technique  by  which  the  IQ  method 
can  handle  transmission  lines  will  be  briefly  discussed. 

Consider  the  lossless  transmission  line  of  Figo  8.5a.   It  was  found 
that  the  m.ost  convenient  equivalent  circuit  is  that  shown  in  Fig,  8o5b.   Z  is 
the  characteristic  im_pedance  of  the  line  and  the  forward  and  reverse  currents 
at  the  sending  and  receiving  end  are  defined  as  follows: 
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Figure  8.4.   Asymmetric  Flipflop  Example, 


-55- 


U  I. 


I 

■o 


>z< 


CHARACTERISTIC      IMP  =  Zq 
VELOCITY    OF  PROPAGATION  =V, 


>Zl 


T  =L 
Q)  ACTUAL  TRANSMISSION  LINE 


b)   EQUIVALENT  CIRCUIT 


[I    NODE    10 


I 


NODES 


l'-©h    a,p' 


?  ? 

I  I 

I  I 

I  I 


[*      NODE  9 


J 


ws 


?  ? 

I  I 

I  I 

I  I 


NODE  5 


? 


? 


NODE  4 


I  i 

I  I 

i  I 


C)  A  TRANSMISSION    LINE  IN  A  CIRCUIT 


Figure   8.5«      Transmission  Line   and  Equivalent 
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3 
1^=1^   at  the  sending  end 


■p 
T  =  T    at  the  receiving  end 
F    F 


■□ 

I  =  I    at  the  receiving  end 


g 
I  -   I    at  the  sending  end 
R    R 


(8.31) 


We  therefore  have  the  following  relations 


13=1^.  I^      and      1,  =  1^^  .  1^  (8.32) 


I^(t)  =  l|(t  -T) 


4(t)  =  l^(t 


(8.33) 


By  imposing  boundary  conditions  so  that  the  reflection  coefficients 
are  correctly  defined  we  get 


S   ^S 

R  '  ^OS 


i;  =  i: .  I 


(8.3^) 
^  =  I?  -  I 

R    "T     OR 


q  "P 

From  the  previous  set  of  definitions  we  see  that  I  and  I^  are  the  delayed 
forward  and  reverse  currents  respectively  and  therefore  are  known  quantities 
at  any  time  instance.   If  the  equivalent  impedances  of  the  line  at  the  sending 
and  receiving  ends  are  treated  as  normal  branches  by  the  IQ  method,  I^g  and  I^^^ 

will  be  known  after  the  current  distribution  cycle.   Consequently  the  currents 

S       R 
I^  and  1^   which  are  to  be  entered  into  the  two  delay  tables  are  easily  calcu- 

F      R 

R      S 
lated.   These  emerge  after  a  time  T  as  the  values  for  I^  and  I^  respectively. 
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The  above  can  also  be  adapted  to  lossy  transmission  lines^  the 
difference  being  that  now  the  data  in  the  delay  tables  must  be  processed  in 
some  manner  Instead  of  being  merely  stored.   Since  this  is  out  of  the  domain  of 
this  report,  the  subject  will  not  be  discussed  further. 

Although  the  above  relations  for  the  transmission  line  are  very 
simple;,  two  complications  enter  during  implement ati on „   Since  the  delay  tables 
must  be  of  fixed  size  there  is  a  problem  of  matching  the  fixed  At  of  the  line 
to  the  variable  time  increment  which  is  possible  in  the  IQ  method.   This  problem 
is  fairly  easily  solved  by  integrating  the  charge  entering  the  line  when  the 
line  At  is  larger  than  the  IQ  time  increment  and  by  "filling"  the  spaces  in 
the  table  which  are  not  used  when  the  line  At  is  smaller  than  the  IQ  time 
increment.   The  second  problem  is  one  of  initializing  the  delay  tables.   Since 
points  1  and  3,  and  2  and  k   in  Fig.  8.5  are  parts  of  the  same  node  at  dc,  they 
should  also  be  treated  by  the  IQ  method  as  the  same  nodes  but  this  presents 
some  problems  which  have  not  yet  been  completely  resolved  and  which  are  directly 
related  to  how  the  main  IQ  method  is  programmed.   However,  this  problem  also 
does  not  seem  to  be  insurmountable. 


9.   REALIZATION  OF  THE  IQ  METHOD 

9.1  The  Program 

To  verify  some  of  the  ideas  contained  in  the  above,  a  program  based 
on  the  IQ  method  was  written.   Although  the  time  element  prevented  the  inclusion 
of  all  of  the  techniques  developed,,  the  basic  ideas  of  sections  2  and  3  were 
all  included.   In  addition  a  method  for  handling  nonlinear  two-terminal  devices 
such  as  tunnel  diodes  and  the  method  for  analyzing  transmission  lines  were 
incorporated  into  the  program.   Since  some  major  changes  in  the  organization  of 
the  program  would  be  required  to  include  efficiently  the  voltage-controlled 
subnode  method^,  the  addition  of  this  and  hence  the  model  for  the  diode  and 
transistor  has  not  yet  been  completed.   The  program  was  written  in  FORTMN  II 
for  the  IBM  709^  system  at  the  Digital  Computer  Laboratory,  University  of 
Illinois , 

The  input  part  of  the  program  is  written  in  such  a  manner  that  the 
sequence  in  which  the  nodes  are  read  automatically  defines  their  level.   For 
each  node  the  numbers  of  all  the  branches  below  the  node  are  also  read  from  the 
tape.   The  branch  data  read  in  by  the  program  consists  of  the  node  number  below 
the  branch,  the  four  constants  of  the  branch  and  if  it  has  a  nonlinear  element, 
the  type  and  number  of  the  nonlinearity.   The  rest  of  the  data  is  read  by  the 
program  in  a  standard  manner. 

The  main  part  of  the  program  consists  of  routines  for  inputting  and 
outputting  data,  incrementing  circuit  variables  and  for  the  actual  calculation 
of  the  response,   A  simplified  flow  chart  of  the  latter  is  shown  in  Fig,  A-1 
of  the  appendix.   After  initializing  of  the  table,  the  dc  initial  conditions 
are  evaluated.   Since  this  uses  the  same  routines  as  the  transient  analysis 
with  minor  changes,  the  program  is  put  in  "Mode  l"  during  the  dc  calculations. 
While  in  this  mode,  the  subroutines  will  replace  all  capacitances  by  open 
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circuits  and  all  inducatances  by  shorts,  which  leads  to  a  relaxed  condition 
when  the  dc  analysis  is  completed,.   Other  initial  conditions  could  be  easily 
specified  by  making  minor  changes  in  the  program,  but  these  have  not  yet  been 
programmed „   Initially  all  sources  have  zero  values;  then  the  program  "turns 
them  on"  by  successively  increasing  the  values  of  the  sources  until  they  have 
their  desired  values,  meanwhile  continually  solving  the  circuit  to  prevent  large 
discontinuities.   This  process  does  not  seem  to  take  excessive  time;  yet  it 
allows  bistable  circuits  to  be  set  as  desired.   This  part  of  the  program  in 
itself  provides  a  means  for  the  dc  analysis  of  nonlinear  circuits. 

When  the  dc  initialization  is  completed,  Mode  2  is  set  and  the 
transient  analysis  is  started.   The  control  of  the  printing,  time  incrementing 
and  impedance  calculation  is  shown  in  the  right  half  of  Fig.  A-l,   QV  is  the 
name  of  the  subroutine  which  distributes  the  charge  and  calculates  the  node 
voltages,  IMP  calculates  the  impedance  factors  and  ALTERT  changes  the  time 
increment  and  checks  if  the  present  set  of  values  is  to  be  printed.   The  sim- 
plified flow  charts  of  IMP  and  ALTERT  are  found  in  Figs,  A-k   and  A-5„   The  terms 
found  in  Fig,  A-5  are  defined  by  the  following^ 

DT   --Present  time  increment 

TMAX  --Time  at  which  analysis  is  stopped 

DTM  --Minimum  time  increment  allowed 

DVM  --Maximum  node  voltage  increase  during  cycle 

DW  --Maximum  node  voltage  increase  allowed 

LVL     --If  DVM  is  less  than  this  the  time  increment  can  be  doubled 

It  is  to  be  emphasized  again  that  the  flow  charts  do  not  represent  the  subrou- 
tines exactly;  their  purpose  is  merely  to  present  some  of  the  major  computational 
steps  in  the  subroutines. 
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Figures  A-2  and  A- 3  contain  the  charge  distribution  and  node  voltage 
calculation  loops  which  will  require  the  most  computer  time  during  the  analysis. 
The  flow  charts  in  this  case  describe  the  steps  discussed  in  section  2  fairly 
completely,  but  again  minor  points  are  not  included.   As  soon  as  QV  is  entered, 
it  calls  a  subroutine  ADJERR  (Fig,  A-6)  which  changes  the  adaptive  factor  v 
discussed  in  section  5„2,   At  present  a  slightly  larger  than  linear  extrapola- 
tion formula  is  used  to  minimize  the  error  in  future  numerical  cycles.   This 
is  done  only  when  there  is  a  definite  increasing  trend  in  the  error  over  the 
last  three  cycles;  otherwise  an  average  is  used  and  if  the  error  is  oscillating 
between  positive  and  negative  values  no  correction  is  made  at  all. 

At  present  the  decrease  of  the  time  increment  is  based  on  how  large 
the  node  voltage  changes  are  rather  than  how  large  the  circuit  error  is.   This 
is  a  result  of  earlier  indecision  about  which  node  voltage  calculation  criterion 
to  use.   Experience  with  the  program  has  indicated  the  time  interval  change 
should  rely  more  on  the  circuit  error, 

9 , 2   Sample  Results 

To  demonstrate  some  of  the  properties  of  the  IQ  method,  four  circuits 
and  their  response  will  be  presented  here.  The  first,  that  of  Fig,  A-7,  illus- 
trates the  effect  of  v  on  the  circuit  error.   The  maximum  allowed  error  was 

a 

set  arbitrarily  high  so  that  the  method  would  not  decrease  the  size  of  the 
increments.   When  v  =  0,  we  see  that  the  error  in  the  initial  two  cycles  after 
the  application  of  the  10~volt  step  is  very  small  but  increases  slowly^  as 
discussed  in  section  5,2,   The  vast  improvement  by  correctly  changing  v^  is 
noted,  confirming  the  discussion  of  that  section. 

In  the  second  example  a  tunnel  diode  monostable  (Fig,  A-8)  is  analyzed, 
The  tunnel  diode  curve  used  is  composed  of  three  cubics  for  the  tunneling, 
valley  and  diffusion  regions  of  the  diode  which  were  found  by  a  least-square 


■61- 


method  of  curve  fittings   The  same  data  was  processed  by  the  IQ  method  and  by 
the  second-order  Runge-Kutta  method,  both  of  which  gave  the  same  response, 
agreeing  within  three  millivolts  at  any  given  time.   The  circuit  error  in  the 
IQ  method  was  always  less  than  one  millivolt  during  the  fast  transients  and 
even  less  at  other  times o 

In  the  third  example  a  tunnel  diode  pulse  amplifier  and  waveform 
reshaper  is  considered^   The  tunnel  diode  on  the  input  side  had  a  peak  current 
of  two  milliamperes  and  that  on  the  output  a  f ive-milliampere  peak.   The  two- 
milliampere  tunnel  diode  triggers  the  five-m.illiamper e  diode  which  then  supplies 
the  output  power „   The  effect  of  the  coupling  resistance  R  is  observed  in  the 
two  response  diagrams „   Since  the  average  circuit  error  during  the  response  was 
less  than  O.5  millivolt,  the  waveforms  can  be  expected  to  be  fairly  accurate „ 

As  a  final  example,  the  reflections  from  an  unbiased  tunnel  diode 
terminating  a  transmission  line  (Fig,  A-IO)  are  analyzed  by  the  IQ  method.   The 
input  pulse  and  the  reflection  appear  in  the  v  versus  t  plot.   The  initial 
negatively-reflected  voltage  occurs  because  of  the  low  tunnel  diode  resistance 
in  the  tunneling  region.   After  the  peak  is  passed  the  resistance  becomes  nega- 
tive, then  very  high,  resulting  in  a  positive  voltage  reflection  coefficient. 
The  final  spike  results  from  the  capacitance  discharging  into  the  line.   Since 
the  sending  end  is  correctly  matched,  it  produces  no  reflections, 

9-3  Efficiency 

The  program  that  is  presently  operative  has  been  written  in  a  manner 
so  as  to  allow  easy  changes;  efficiency  was  not  a  primary  objective  and  con- 
sequently not  much  information  can  be  gained  by  comparing  the  speed  of  computa- 
tion of  the  present  program  with  that  of  other  methods,   A  more  significant 
approach,  an  examination  of  the  relative  time  spent  in  different  parts  of  the 
program,  indicates  that  the  calculation  of  the  branch  currents  and  voltages  in 
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QV  by  far  uses  the  most  time.   This  then  leads  to  the  conclusion  that  the  time 
required  to  analyze  a  circuit  will  "be  proportional  to  the  number  of  branches 
it  haso   To  get  an  order  of  magnitude  estimate  of  the  speed  of  the  IQ  method 
on  a  computer  we  need  therefore  only  consider  the  operations  required  in  the 
branch  loops  of  the  program. 

For  ILLIAC  11^  for  instance,  it  was  found  that  about  20  additions^ 
30  multiplications  and  50  access  times  will  be  required  per  branch  calculation. 
Based  on  times  of  3-5  [J^sec,    6  (asec^,  and  1,8  lasec  respectively  for  these  three^ 
this  comes  to  a  total  of  350  |-isec  per  branch.   Since  branches  with  nonlinear 
elements  will  require  more  time^  we  can  roughJ-y  say  that  the  average  time  per 
branch  calculation  is  500  lasec.   Therefore  a  20-branch  circuit  requiring  1000 
steps  for  the  solution  should  take  about  10  seconds. 

Two  possible  methods  of  increasing  the  speed  of  the  IQ  method  should 
be  mentioned.   Firsts  in  most  circuits  there  are  many  branches  which  are  con- 
nected to  ground;,  the  reference  node  for  the  IQ  method.   Since  these  often 
contain  only  single  elements  such  as  resistances  or  capacitance S;,  more  simplified 
versions  of  the  correction  formula^  equation  2, l6^  could  be  used  for  these  special 
cases.   Since  sometimes  as  many  as  half  of  the  branches  will  be  connected  to  the 
reference  node^  this  represents  a  significant  increase  in  speed. 

Second,  for  large  circuits  it  is  possible  that  one  end  is  undergoing 
transients  but  the  other  is  still  in  a  steady  state.   Therefore  the  two  parts 
of  the  circuit  can  be  analyzed  with  different  time  increments.   Although  this 
will  lead  to  problems  in  matching  the  time  increments  the  result  should  be 
especially  rewarding  during  the  analysis  of  large  transistor  switching  circuits. 


10 „   CONCLUSIONS 

At  present  not  enough  experience  has  been  had  with  the  IQ  method  to 
draw  very  definite  conclusions  about  it.   But;,  to  summarize,  we  can  restate 
some  of  its  advantages  and  disadvantages „   Among  the  advantages  are  the 
following: 

1)  An  absolute  error  is  available  which  the  IQ  method  uses 
to  make  corrections „ 

2)  An  error  correction  mechanism  is  provided  by  the  method 
for  nonlinear  elements,  which  eliminates  the  need  for 
repetitive  relaxation  formulas o 

3)  The  method  is  particularly  suited  for  analyzing  diode 

and  transistor  circuits  since  the  provided  error  correction 
mechanism  eliminates  the  need  to  integrate  the  transistor 
and  diode  equations  separately  and  match  the  solutions 
with  the  rest  of  the  circuit  periodlcallyo   A  second 
merit  Is  that  no  further  assumptions,  aside  from  those 
used  to  get  the  Ebers  and  Moll  equations,  need  to  be  made 
because  of  the  convenient  manner  in  which  the  relevant 
transistor  quantities  can  be  calculated. 
h)      The  efficiency  seems  reasonable  and  the  time  required  to 
analyze  a  circuit  is  linearly  dependent  on  the  number  of 
branches  in  the  circuit „ 

5)  Programming  is  not  too  complicated  since  no  matrices  must 
be  prepared  and  the  initial  conditions  and  transients  can 
be  calculated  by  the  same  routine. 

6)  The  IQ  method  is  very  much  directed  toward  giving  the 
user  insight  into  his  circuit.   The  method  deals  with 
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quantities  he  can  visualize  and  allows  him  to  control 
more  directly  how  the  program  handles  his  circuit. 

The  chief  disadvantage  is  that  the  method  will  not  handle  all  circuits 
with  equal  ease  since  it  was  developed  for  a  circuit  with  a  restricted  topology. 
This  will  require  some  work  from  the  user  before  he  has  his  circuit  in  a  form 
in  which  the  method  will  handle  it  most  efficiently.   How  serious  these  limita- 
tions can  be  is  not  yet  know. 

It  is  felt  that  at  present  the  IQ  method  is  still  in  too  much  of  an 
untested  state  to  allow  a  detailed  comparison  with  other  numerical  methods  for 
circuit  analysis.   But,  considering  the  above  advantages  of  the  IQ,  method^  it 
is  also  felt  that  it  will  certainly  prove  to  be  a  useful  tool  and  therefore 
deserves  further  analysis  and  refinement. 

As  a  final  note,  it  should  be  pointed  out  that  the  IQ  method  developed 
in  this  report  is  a  consequence  of  the  chosen  topology  and  set  of  differences. 
It  is  possible  that  other  incremental  charge  methods  with  better  characteristics 
can  be  formulated  by  applying  the  same  basic  principles  and  deriving  parallel 
formulas , 
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Figure  A-lo   Flow  Chart  for  Main  Routine, 
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Figure  A-2 ,      Flow   Chart  for  the   Charge   Distribution  Cycle. 
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Figure  A-3.      Flow  Chart  for  Node  Voltage   Calculation. 
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Figure  A-4.      Flow   Chart   for   Impedance  Factor  Calculations 
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Figure  A-5.      Flow  Chart   for  At  Alteration  Subroutine 
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Figure  A-8o   Tunnel  Diode  Monostable  Response  by  the  IQ,  Method, 
and  by  the  Second-Order  Runge-Kutta  Method, 
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Figure  A-IG,      Reflections   from  an  Unbiased  Tunnel  Diode „ 


